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A theoretical description of the radial density profile for charged particles with Yukawa inter- 
action in a harmonic trap is described. At strong Coulomb coupling shell structure is observed 
in both computer simulations and experiments. Correlations responsible for such shell structure 
are described here using a recently developed model based in density functional theory. A wide 
range of particle number, Coulomb coupling, and screening lengths is considered within the fluid 
phase. A hypernetted chain approximation shows the formation of shell structure, but fails to give 
quantitative agreement with Monte Carlo simulation results at strong coupling. Significantly bet- 
ter agreement is obtained within the hypernetted chain structure using a renormalized coupling 
constant, representing bridge function corrections. 
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I. INTRODUCTION 

Spatially confined charged particles have attracted 
growing interest. Examples include electrons in quan- 
tum dots [1 , ions in Penning and Paul traps [21 13] and 
the mesoscopic charges of dusty plasmas [il [5] . In partic- 
ular, three-dimensional classical spherical plasmas have 
been produced in ion systems [6j and more recently in 
dusty plasmas [7]. The structural and dynamic proper- 
ties of these systems continue to attract the interest of 
many groups in various fields; e.g., [SHID]- 

At sufficiently strong coupling these systems form con- 
centric shells which are well reproduced by simulations, 
c.f. [TTUI3] and references therein. The objective here is 
to provide a theoretical analysis to complement these re- 
sults from simulation and experiments, for a better phys- 
ical understanding of the underlying mechanisms. For 
harmonically confined particles with Coulomb interaction 
such a theory of shell formation as a function of tempera- 
ture (inverse coupling) was derived recently using classi- 
cal density functional theory (DFT) p4l[T6] . However, a 
special property of dusty plasmas is the screening of the 
pair interaction. The theoretical description is extended 
here to describe spherically trapped strongly correlated 
particles with Yukawa interaction for such dusty plasmas. 

The state conditions are specified by three dimension- 
less parameters: particle number N , coupling constant F 
(defined below), and dimensionless screening parameter 
K*. The ranges of values considered are 15 < A'^ < 500, 
< F < 100, and < k* < 1. The primary focus here 
is on shell formation as a function of these parameters. 
Only the equilibrium fluid phase is considered (rotational 
invariance) so that shell structure is reflected in the ra- 



dial dependence of the density of confined charges. The 
average density is defined as a multi-dimensional config- 
uration integral in the canonical ensemble, which can be 
evaluated by Monte Carlo simulation. New simulations 
are provided here as a means to determine the accuracy of 
theoretical approximations. The system, dimensionless 
units, and adaptation of the hypernetted chain (HNC) 
theory introduced in 'T^ for confined Coulomb charges 
are described in Section |II A| The density profiles from 
the HNC approximation are compared to simulations in 
Section |IH[ It is found that the formation of shells, as 
well as their location and populations, is well described 
by the HNC approximation, but the shell maxima and 
widths show large discrepancies for F > 10 and the er- 
rors increase with increasing k*. The primary qualita- 
tive effect of screening is to shift the shells toward the 
center and decrease the overall volume. An "adjusted" 
hypernetted chain approximation (AHNC) is considered 
in Section |IV[ This is based on a model for the bridge 
function corrections to HNC first introduced by Ng [17] 
for the pair distribution function of a Coulomb one com- 
ponent plasma (OCP). It has the property of preserving 
the form of the HNC equations with only a renormaliza- 
tion of F to some larger effective value. It is shown that 
the same method applies to the Yukawa OCP for accu- 
rate pair correlations even at very strong coupling, and 
the approach is then applied to the bridge corrections to 
the equation for the density profile. An optimized renor- 
malization for the density leads to excellent results for 
the Coulomb case (e.g., F < 100, < 500). For Yukawa 
systems this approach is very useful as well for the radial 
distribution function, but it is somewhat more limited 
for the density profile at larger values of k* and N . 
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II. THEORY AND SIMULATION 

A. Model and units 

The system is comprised of N identical charges inter- 
acting pairwise via a Yukawa potential, confined by a 
harmonic potential centered at the origin. The Hamilto- 
nian is 



H = 



El 



2 2 2 

-mVj + -muj r,- 
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Here m is the mass, uj is the angular frequency measuring 
the strength of the confinement, and Y^,^fi are the position 
and velocity of charge i. The Yukawa interaction is 



(2) 



where rij = Ir^ — rj|, q is the particle charge, and n 
is an inverse screening length. The physical origin of 
this screening length is described elsewhere [S] and will 
not be discussed here. The primary property of inter- 
est here is the local density of charges in the trap at 
equilibrium. For the classical canonical ensemble its di- 
mensionless form is given by 

f rlr* rlr* p^^*('^1'---''^*n) 
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Here, r* = r^/ro, /3 = l/kgT is the inverse temperature, 
r = I3q^ /r^j is the Coulomb coupling constant, and k* = 
Krp. The usual choice for the length scale tq is the ion 
sphere radius, or mean distance between charges, given 

by 



47rrj^n/3 = 1, 



(5) 



where n is a characteristic spatially averaged density to 
be chosen for convenience. Here it is chosen to simplify 
the Hamiltonian by the condition 



2^2 = 2' 



3mcj2 



(6) 



This is not the average density for the Yukawa particles 
in the trap ut — N/Vr, where the volume Vr — 4:TtRj./3 
is defined by the maximum radius Rt at which the force 
on a charge due to the trap is equal to that of all other 
charges 

mw^i?T = / dr' 2 (1 + ^ - r'l) nrir'). 

J IRt-i-'I 

(7) 



It follows that ut = n for k = and Ht > n for k > 0. 
The solution to ([t]) is discussed further below. 

The dimensionless trap potential energy is now a func- 
tion of two dimensionless parameters, the Coulomb cou- 
pling strength T and the screening parameter k* 
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and the dimensionless density profile n*(r*) can be ob- 
tained numerically from a Metropolis Monte Carlo sim- 
ulation for given F, n* , and N . 



B. Theory and approximations 

A formal representation of the average density profile 
was developed within density functional theory in refer- 
ence |15j . That analysis applies here as well, with only 
the replacement of the Coulomb potential by the Yukawa 
potential. First, the density is represented in terms of a 
dimensionless effective potential U* (r*) 



n*(r*) EE N 



-ru*{r') 
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Here N denotes the average number of particles in the 
trap, since the theory is formulated in the grand canoni- 
cal ensemble. The effective potential obeys the equation 



r*'|) 



+B{r*). 
(10) 

The function c(r*) is proportional to the direct correla- 
tion function for a uniform one component plasma (OCP) 
of Yukawa charges 



(11) 



which is related to the OCP radial distribution function 
5ocp(^*) by the Ornstcin-Zernike equation 

<7ocp(r-*)-l = cocp(»-*)+ (12) 
+ fT^cp / c^r*' [gocp(r*') - 1] cocp(|r* - r*'|). 



(13) 



Finally, gocp (''* ) is determined from the equation 



ln5ocp(?'*) 



n, 



OCP j dr*' [gocp(r*') - 1] cqcp (|r* - r*'|) 
rSocp (r*). 



The second term of ( 10 ) describes the effect of cor- 



relations among particles in the trap in terms of the 
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corresponding correlations among particles in the uni- 
form OCP. The last term B [r*) corrects this approxi- 
mate treatment of correlations and is known as a bridge 
function. Similarly, i?ocp (r*) is the bridge function for 
<?ocp('"*) [H]- To optimize this contribution of OCP cor- 
relations, the density of the trap is matched to that of 
the OCP 



'OCP 



(14) 



For given TV the trap density is fixed by the volume of the 
trap, whose radius Rt must be calculated from ([t]). An 
approximate evaluation for the ground state has been 
discussed elsewhere [T^, with the result that it is the 
unique positive, real solution to 



(1 + K*R*){N - 1) + R*^ + K*R 
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In all of the following, is determined in this way for 
each K* . 

The above equations ^ - (13) are still exact, but re- 
quire specification of the bridge functions. The simplest 
approximation is the neglect of the bridge functions, lead- 
ing to the hypernetted chain approximation (HNC) 



/ dr*'e~'"'^SNc ('"*') 
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(16) 
(17) 
), 

(18) 



gii-Ncir*) - 1 = CHNC(''*) + 

+n*^ J dr*'(gHNc(r*')-l)cHNc(|r*-r*'|). 



This is a closed set of equations for C^hnC' 5hnc(^*)i and 
Chnc- Note that the determination of 5hnc(''*) and chnc 
is independent of the trap density calculation. 

It is well known that the HNC approximation for 
the OCP properties is a good approximation except for 
strong coupling where the bridge function i?ocp becomes 
important. However, the results below for the trap den- 
sity show that the trap bridge function can be important 
even at moderate coupling. It is therefore necessary to go 
beyond HNC and find an approximation for the bridge 
functions. This is described below. 



III. RESULTS: HNC APPROXIMATION 

Correlations in the HNC approximation are described 
by Chnc- For weak coupling, P < 1, chnc ~^ e"** ^ 1^* ■ 
This is the "mean field" limit. Fig. [TJshows a comparison 
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FIG. 1: Mean field results for the density profile (lines) for a 
Coulomb and Yukawa OCP in a spherical trap for (a) P = 1 
and (b) F = 5 for A'' = 100 compared with Monte Carlo 
simulations (symbols). 



of this mean field description with Monte Carlo simula- 
tion results at moderate coupling, F = 1 and 5 for several 
values of k* . As might be expected, there is reasonable 
agreement at F = 1, but emergence of an outer shell is 
evident at F = 5, which cannot be reproduced by the 
mean field theory. Evidently, here it is necessary to cal- 
culate the correlations of chnc through the full coupled 
set of equations ( [l6|) - (18). 

Figures [2^) and [2]3) show chnc as a function of r* 
for K* — Q and 1. In both cases the deviation from the 
mean field limit increases for stronger coupling, creating 
a "correlation hole" for r* < 1. The effects of these 
correlations on the trap density profile are illustrated for 
several values of the screening parameter k* in Fig. [3] 
at F = 50, N ~ 100. It is seen that increased screening 
tends to compress the system [TT] and enhance the shell 
structure. 

The quality of HNC is tested by comparison to Monte 
Carlo simulations. This is illustrated for = 100 and 
F = 10, 40, 100 with k* = in Fig. |4]a) and with k* = 1 
in Fig. |4]d). HNC is a poor approximation at r* = 
which results in overall poor results for small particle 
numbers A^ < 10. This error appears periodically with 
the creation of each new shell and is small if no particle 
is at the center. For k* = the shell locations match 
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FIG. 2; Various direct correlation functions for (a) Coulomb 
interaction and (b) Yukawa interaction with k* = 1. The top 
curve is the mean field value. 
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FIG. 4: Density profile for A'^ = 100 particles and various F 
values: comparison of HNC results (solid lines) with Monte 
Carlo (symbols) for (a) Coulomb and (b) Yukawa interaction 
with K* = 1. 



require going beyond the HNC approximation. 
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FIG. 3: HNC density profile for a Yukawa system with various 
K* at r = 50 and iV = 100. 



IV. ADJUSTED HNC 

In a recent analysis for Coulomb systems in a spheri- 
cal trap it was also observed that the HNC approxima- 
tion gives the correct location and population of shells 
[HI [TB] , which depend only weakly on F. For Yukawa 
systems, these properties become less accurate with in- 
creasing At*. For both Coulomb and Yukawa, the ampli- 
tude and width depend strongly on F and are underes- 
timated by the HNC approximation at strong coupling. 
This suggests that increasing the coupling constant alone 
would increase the accuracy of HNC. 



A. Pair distribution function 



well the simulation data, while increasing k* leads to 
decreased accuracy for the inner shells. The effect is small 
up to K* = 0.5. Figure [5] compares HNC and Monte Carlo 
results with N — 15, 125, and 500 for Coulomb charges 
in Fig. [5^), and for Yukawa charges with k* = 1 in Fig. 
[5}d). The shell populations (not shown) and locations are 
nearly independent of F, as seen in MC simulations [2Ql 
m], and are given accurately by HNC. However, the peak 
heights and widths for the shells are poorly predicted and 



This failure of HNC for strong coupling has been stud- 
ied in some detail for the calculation of the Coulomb 
5ocp(^)- Among the earliest investigations is that of Ng 
[T7] who observed that the HNC peak positions are given 
accurately for strong coupling, but not the amplitudes 
and widths. He corrected the HNC by representing the 
bridge function of ( 13 1 in the form 
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FIG. 5: Density profile for F — 20: comparison of HNC results 
(lines) with Monte Carlo results (symbols) for (a) Coulomb 
and (b) Yukawa interaction with n* — 1. 
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FIG. 6: Comparison of AHNC results for the pair distribution 
function (/(r*) with simulations for (a) Coulomb and (b) k* = 
2. The values of F' in (a) are 12.5, 57 and 160, in (b) 10, 130 
and 480. 



where A(r) is a chosen function of F and V{r*^ is the 
Coulomb potential. This particular choice was not ob- 
tained from any theoretical analysis, but rather because 
it leads back to the HNC form with a renormaUzed cou- 
pling constant F' = [1 -I- A(F)]F. This approach will be 
referred to as the adjusted HNC (AHNC). It was shown 
that an accurate prediction of gocp(^) could be obtained 
over the entire fluid domain with the choice 



A(F) ^ Ajvg(r) = 0.6erf(0.024r). 



(20) 



Subsequent theoretical studies of the Coulomb bridge 
function by Rosenfeld and Ashcroft [22^ , indicated that it 
has a " universal" form and hence could be represented by 
the corresponding hard sphere bridge function for which 
an accurate parametrization is known. Although consid- 
erably more complex to implement computationally, it 
also gives a very good representation for goc^ir)- Fur- 
thermore, it has an important thermodynamic consis- 
tency not shared by the HNC or AHNC approximations. 
Evidently, the functional form ( 19 1 represents the actual 



bridge function for the relevant range of r needed to de- 
termine gocp (j" ) (the numerical difficulty of determining 
Bqcp (r*) precisely from 5ocp(^*) is discussed by Poll et 
al. [53]) Due to its simplicity and the direct interpreta- 
tion as a renormalization of the coupling strength the 
AHNC will be used here as the means to improve the 
HNC approximation. 

It remains to show how the bridge function should be 
chosen for the Yukawa potential. An empirical choice has 
been suggested in the form [24] 



B 



V(r^) -> Socp 



(r*)e- 



V2 



(21) 



where Bqcp (r* ) is the Coulomb bridge function. This 
gives very good results for (7y(^*) when Bqcp (r*) is 
approximated by the corresponding hard sphere bridge 
function, as suggested by Rosenfeld and Ashcroft. In 
contrast the AHNC for the Yukawa potential is obtained 
from pgj) 



SY(r*)^A(F)/?- 



(22) 

where A(F) is the same form as 120] ) for the OCP 

A(F) ^ XNgic{K*)r) = 0.6erf(c(K*)r). (23) 

The constant c{k*) is adjusted for each k* , with the Ng 
value c(0) = 0.024. This Yukawa AHNC leads to the 
HNC form ( [l7| , but with a renormalized coupling con- 
stant 

lngAHNc(?'*) = -r' h 



+ n*^J dr*'[5AHNc(r-*')-l]cAHNc(|r*-r*'|)(24) 



Figures [6p,) and |6p) show the excellent agreement be- 
tween AHNC and molecular dynamics even at very 
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strong coupling for both k* = and k* = 2. (Note that 
r and K* used here refer to a length unit tq defined by 
Eq. ([5]) with the OCP density). It is interesting to note 
that results for recent MD simulations for different values 
of r and K* can be collapsed in terms of a single effective 
coupling constant F* = T* (F, k*) |25j . In summary the 
AHNC for g{r*) proposed by Ng for the Coulomb OCP 
works as well for the strongly coupled Yukawa plasma. 



B. Density profile 

With the results for the homogeneous OCP pair dis- 
tribution function as a guide, a similar representation 
is considered for the bridge function B{r\n) of the trap 
density profile ( 10 ) 



Bir) A(r)/3Vb(r*). 



(25) 



Here [3Vo is the trap potential rr*^/2, restoring the HNC 
form ^ and ^ with a renormalized F'. An initial ap- 
proach would be to use the same renormalization function 
A(F) as obtained in the optimization of (Jahnc- This im- 
proves the accuracy for coupling constants up to F ~ 40. 
To include stronger coupling it is necessary to choose 
a different renormalization function A(F) when calculat- 
ing the trap density profile. Although equations ( 19 ) 



and ( 25 ) formally allow for separate specifications of the 



renormalization function A(r) for the OCP and trap sys- 
tems, results show that the same function A(F) must be 
used to determine the trap density profile to agree with 
simulations. That is, A(F) can be determined by fitting 
either the density profile or the pair correlation function, 
but not both. 

The explanation as to why these two approaches (de- 
termining A(r) separately for the systems and as a com- 
mon function) give very different results when using ( 19 ) 
and (251, lies in the relationship between direct correla- 
tion functions at different F. The scaled direct correla- 
tion function (111 is independent of F for F > 10. That 
is, if separate coupling constants for the trap (Ftrap) and 
OCP (Focp) systems were fitted, the direct correlation 
functions are still related by 



c(r; Ftrap) _ c(r;Focp) 



Ptra 



OCP 



(26) 



By considering again ( 10 1, this shows that the two proce- 



dures of using one common coupling constant, and using 
separate coupling constants, are related by scaling the 
number of particles in the trap. As the procedure of 
using a common renormalization function has shown to 
give good results, an equivalent approach involving sep- 
arate renormalization functions is to have the effective 
number of particles in the trap also be dependent on the 
coupling constant so as to correct the discrepancy. It is 
important to note that with this alternative approach, 
although both the coupling constant and particle num- 
ber are scaled, there still is only one fitting parameter 
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FIG. 7: Depend ence of the AHNC parameter A on F and 
K* calculated by (271 at fixed particle number N = 100. The 
renormalized coupling parameter is given by F' — [1 + A(r)]r. 
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FIG. 8: Dependence of the AHNC parameter A on A for the 
Coulomb case. The height of the outermost peak is used to 
obtain A. The renormalized coupling parameter is given by 

r' = [1 + A(r)]r. 



for the trap (which fixes both a scaled coupling constant 
and particle number), and one fitting parameter for the 
OCP. However, the interpretation of the scaled particle 
number for the trap is not clear, as the shell structure 
depends critically on N. 

Therefore we proceed by choosing to fit only the trap 
density profile, and using a common renormalization 
function for both the trap and the OCP systems. 

An appropriate value of A(F) to optimize the density 
profile is obtained by minimizing the square difference 
of the Monte Carlo data rtMci^) and the HNC profile 
nHNc{T\^') with respect to F' 



F' : 



drr^ [nMcir\T) - nHNc('"|r') 



(27) 



Since at small r the HNC density profile is not accurate, 
the difference in the peaks is weighted to the particle 
number in each shell. In practice this effectively fits the 
height and width of the outermost peak. The dependence 
of XiT,N) on F and N is shown in Fig. [t] and Fig. [s] 
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The r dependence for different k* cannot be collapsed to 



a single curve by rescaling F as in ( 23 ) , as the asymptotic 
large k* limits of A(r, N) are now different. 

The quality of the AHNC approximation is again es- 
tablished by comparison to Monte Carlo data, cf. Fig. 
[9j AHNC describes accurately the density profile of 
Coulomb charges for the full rang e of F [Fig. H (a)] and 



10 (a)], while keeping the sim- 



particle numbers N [Fig. 
pie form of HNC. A similar improvement in accuracy 
is observed for the Yukawa system with k* < 0.5 and 
N < 100 [Fig. 9] (b)]. For larger k* errors in the inner 
shells occur and increase with increasing k* and N [Fig. 
[9] (c) and Fig. lO] (b)]. For large particle numbers fit- 
ting A becomes more subjective, depending upon which 
criteria are imposed, e.g. best outermost peak height or 
inner shell heights. Similarly, increasing the renornialized 
coupling constant beyond a certain value is only trading 
agreement from inner to outer shells. 

It is curious that the AHNC procedure works so well 
for pair correlations, both Coulomb and Yukawa, and 
for the Coulomb trap density profile, but fails for the 
Yukawa density profile at large k* and N. One possible 
explanation is the following. The density equation of the 
HNC entails an additional approximation not contained 
in that for the pair correlations, namely that the pair cor- 
relations in the trap can be represented by those of the 
OCP. This can be justified for Coulomb interactions, but 
that argument does not extend to Yukawa interactions. 
At large k* this approximation may no longer hold. In 
addition, the shell structure is enhanced at large k*, and 
the number of shells increases with N. Hence there are 
increased demands on the AHNC to represent more com- 
plex structure. 

There is a qualitative difference between the Coulomb 
and Yukawa cases at large N. In the former case, the 
harmonic trap is exactly equal to the effect of a uniform 
neutralizing background, and the system approaches the 
Coulomb OCP for large N except at the boundaries. 
However, for the Yukawa case the relationship of the trap 
to the neutralizing background no longer holds. Totsuji 
et al. have derived the corresponding confinement poten- 
tial for a Yukawa system [57] . 



V. DISCUSSION 
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FIG. 9: Density profile for A'^ = 100 particles: Comparison of 
AHNC results (lines) and MC data (symbols) for (a) Coulomb 
interaction, and Yukawa interaction with (b) k* — 0.5 and (c) 
K* = 1. The renormalized F' values used are shown in Fig. u\ 



A theoretical description is developed for the shell 
structure of spherically confined Yukawa plasmas. While 
the precise shell occupations are well known from com- 
puter simulations, both for trapped Coulomb, e.g. [121 
[T5] and Yukawa plasmas e.g. [TT], it is desirable to have 
an analytical theory that correctly reproduces these re- 
sults and provides physical insight into the correlation 
properties. Classical density functional theory is the 
proper starting point for this. In particular, it has been 
shown that the HNC approximation is able to provide the 
density profile (the formation, shape, location, and popu- 
lation of shells) accurately for weak to moderate coupling 



(F < 10). However, HNC fails to reproduce the correct 
width of the shells. 

A simple representation of the bridge functions B (cor- 
rections to HNC) called adjusted HNC (AHNC) is able 
to provide quantitative agreement in the case of Coulomb 
interactions for F < 100 and N < 500, indicating that a 
simple renormalization of the HNC is sufficient to capture 
the structural effects of confinement. A similar adjusted 
HNC provides substantial improvement for the isotropi- 
cally trapped Yukawa system as well. While it correctly 
reproduces the shape of the outermost shell(s) that host 
the majority of particles, it is less accurate for the inner 
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shells, in particular with increasing screening parameter 
and particle number. 
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